/*
 * RandomWalks.cpp
 *
 *  Created on: Jan 21, 2010
 *      Author: snirgaz
 */

#include "RandomWalks.h"



using namespace std;
// Constructor

//Destructor
RandomWalks::~RandomWalks() {
}

void RandomWalks::generateLevelPath(MoveParams *moveParams) {
	IntList::const_iterator itd;
	BeadVec::iterator itb;
	int numOfBeads = moveParams->deltas.size(), b;
	ParticlePos jump;
	double sigma, newPos;
	for (itd = moveParams->deltas.begin(), itb = moveParams->currentBeads.begin();
			itb != moveParams->currentBeads.end(); itd++, itb++) {
		// Draw RW jump
		generate(jump->begin(), jump->end(),
					[&]()->double {return normal_.draw();});
		//QmcRnd::normal_.getMulti(QmcParams::dim, jump.getData());
		sigma = preCalcedSigmas_[*itd];
		jump*=sigma;
#ifdef ROT
		jump.rotate();
#endif
		moveParams->beads[*itb].beadPos_+=jump;
		//itb->beadPos_.addBeadPos(jump);
		assert(moveParams->beads[*itb].beadPos.inBoxLimits());
	}
}

void RandomWalks::init(SimulationParameters& qmcParameters, SimGen* rngGen) {
	normal_.init(rngGen, NormalDistParams(0, qmcParameters.sigma.getVal()));
	int j;
	int M=qmcParameters.M.getVal();
	preCalcedSigmas_.resize(M + 2);
	for (int i = 1; i < (M + 1); i++) {
		double total = double(i), alpha = double(i / 2) / total;
		preCalcedSigmas_[i] = sqrt(total * alpha * (1 - alpha));
	}
}

void RandomWalks::generateJumpBeadPos(BeadVec::const_iterator const &beadStart,
		BeadVec::const_iterator const &beadEnd, BeadVec::iterator &beadTau) {
	int deltaRight, deltaLeft, numOfBeads;
	numOfBeads = beadEnd->RWSlice_ - beadStart->RWSlice_ + 1;
	deltaLeft = beadTau->RWSlice_;
	deltaRight = numOfBeads - (deltaLeft + 1);
	generate(beadTau->beadPos_->begin(), beadTau->beadPos_->end(),
			[&]()->double {return normal_.draw();});
	generate(beadTau->jumpPos_->begin(), beadTau->jumpPos_->end(),
				[&]()->double {return normal_.draw();});
	beadTau->beadPos_*=sqrt(double(deltaLeft));
	beadTau->jumpPos_*=sqrt(double(deltaRight));
//	beadTau->beadPos_.mulBeadPos(sqrt(double(deltaLeft)));
//		beadTau->jumpPos_.mulBeadPos(sqrt(double(deltaRight)));
#ifdef ROT
	beadTau->beadPos_.rotate();
	beadTau->jumpPos_.rotate();
#endif
	beadTau->beadPos_+=beadStart->beadPos_;
	beadTau->jumpPos_+=beadEnd->beadPos_;
//	beadTau->beadPos_.addBeadPos(beadStart.beadPos_);
//	beadTau->jumpPos_.addBeadPos(beadEnd.beadPos_);
	assert(beadTau->beadPos_.inBoxLimits());
	assert(beadTau->jumpPos_.inBoxLimits());
}

